suppressWarnings(library(bitops))
suppressWarnings(library(splines))
suppressWarnings(library(gtools))
suppressWarnings(library(gdata))
suppressWarnings(library(caTools))
suppressWarnings(library(grid))
suppressWarnings(library(survival))
suppressWarnings(library(gplots))
suppressWarnings(library(mvtnorm))
suppressWarnings(library(multcomp))
suppressWarnings(library(R2HTML))
suppressWarnings(library(lattice))
suppressWarnings(library(nlme))
suppressWarnings(library(contrast))
#suppressWarnings(library(lme4))

# retrieve args
Args <- commandArgs(TRUE)

#Read in data
statdata <- read.csv(Args[3], header=TRUE, sep=",")

#Copy Args
#model <- Args[4]
model <- as.formula(Args[4])
timeFactor <- Args[5]
subjectFactor <- Args[6]
covariateModel <- Args[7]
covariance <- Args[8]
responseTransform <- Args[9]
covariateTransform <- Args[10]
blockFactors <- Args[11]
showANOVA <- Args[12]
showPRPlot <- Args[13]
showNormPlot <- Args[14]
sig <- 1 - as.numeric(Args[15])
showLSMeans <- Args[16]
dimfact<-length(unique(eval(parse(text = paste("statdata$", timeFactor)))))
primFactor<- "NULL"


covar <- "NULL"
covar2 <- "NULL"
if(covariateModel != "NULL")
{
	covar<-unlist(strsplit(covariateModel, "~"))[1]
	covar2<-unlist(strsplit(covariateModel, "~"))[2]}
#covar

#Setup the html file and associated css file
htmlFile <- sub(".csv", ".html", Args[3]); #determine the file name of the html file
.HTML.file = htmlFile
cssFile <- "r2html.css"
cssFile <- paste("'",cssFile,"'", sep="") #need to enclose in quotes when path has spaces in it
HTMLCSS(CSSfile = cssFile)

#Output HTML header
if (dimfact ==2)
	{
	HTML.title("<bf>SilveR Within-Subject Parametric Analysis (Paired t-test)", HR=1, align="left")
	} else {
	HTML.title("<bf>SilveR Within-Subject Parametric Analysis", HR=1, align="left")
		}
# Setting up the parameters
resp <- unlist(strsplit(Args[4],"~"))[1] #get the response variable from the main model
statdata$subjectzzzzzz<-as.factor(eval(parse(text = paste("statdata$", subjectFactor))))
statdata$Timezzz<-as.factor(eval(parse(text = paste("statdata$", timeFactor))))
statdata<-statdata[order(statdata$subjectzzzzzz, statdata$Timezzz), ]









#Removing illegal charaters

YAxisTitle <-resp
CPXAxisTitle <-timeFactor

if(covariateModel != "NULL")
{
	XAxisTitle<-covar2
}

for (i in 1:10)
{


# Additional characters included Aug 2010 (STB)
YAxisTitle<-sub("ivs_tilde_ivs"	,"~", YAxisTitle) 
YAxisTitle<-sub("ivs_star_ivs"	,"*", YAxisTitle) 
YAxisTitle<-sub("ivs_plus_ivs"	,"+", YAxisTitle) 

YAxisTitle<-sub("ivs_sp_ivs"	," ", YAxisTitle) 
YAxisTitle<-sub("ivs_ob_ivs"	,"(", YAxisTitle) 
YAxisTitle<-sub("ivs_cb_ivs"	,")", YAxisTitle) 
YAxisTitle<-sub("ivs_div_ivs"	,"/", YAxisTitle) 
YAxisTitle<-sub("ivs_pc_ivs"	,"%", YAxisTitle) 
YAxisTitle<-sub("ivs_hash_ivs"	,"#", YAxisTitle) 
YAxisTitle<-sub("ivs_pt_ivs"	,".", YAxisTitle) 
YAxisTitle<-sub("ivs_hyphen_ivs","-", YAxisTitle) 
YAxisTitle<-sub("ivs_at_ivs"	,"@", YAxisTitle) 
YAxisTitle<-sub("ivs_colon_ivs"	,":", YAxisTitle) 
YAxisTitle<-sub("ivs_exclam_ivs","!", YAxisTitle) 
YAxisTitle<-sub("ivs_quote_ivs"	,"`", YAxisTitle) 
YAxisTitle<-sub("ivs_pound_ivs"	,"£", YAxisTitle) 
YAxisTitle<-sub("ivs_dollar_ivs","$", YAxisTitle) 
YAxisTitle<-sub("ivs_hat_ivs"	,"^", YAxisTitle) 
YAxisTitle<-sub("ivs_amper_ivs"	,"&", YAxisTitle) 
YAxisTitle<-sub("ivs_obrace_ivs","{", YAxisTitle) 
YAxisTitle<-sub("ivs_cbrace_ivs","}", YAxisTitle) 
YAxisTitle<-sub("ivs_semi_ivs"	,";", YAxisTitle) 
YAxisTitle<-sub("ivs_pipe_ivs"	,"|", YAxisTitle) 
YAxisTitle<-sub("ivs_slash_ivs"	,"\\", YAxisTitle) 
YAxisTitle<-sub("ivs_osb_ivs"	,"[", YAxisTitle) 
YAxisTitle<-sub("ivs_csb_ivs"	,"]", YAxisTitle) 
YAxisTitle<-sub("ivs_eq_ivs"	,"=", YAxisTitle) 
YAxisTitle<-sub("ivs_lt_ivs"	,"<", YAxisTitle) 
YAxisTitle<-sub("ivs_gt_ivs"	,">", YAxisTitle) 
YAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", YAxisTitle) 



# Additional characters included Aug 2010 (STB)
CPXAxisTitle<-sub("ivs_tilde_ivs"	,"~", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_star_ivs"	,"*", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_plus_ivs"	,"+", CPXAxisTitle) 

CPXAxisTitle<-sub("ivs_sp_ivs"		," ", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_ob_ivs"		,"(", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_cb_ivs"		,")", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_div_ivs"		,"/", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_pc_ivs"		,"%", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_hash_ivs"	,"#", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_pt_ivs"		,".", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_hyphen_ivs"	,"-", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_at_ivs"		,"@", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_colon_ivs"	,":", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_exclam_ivs"	,"!", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_quote_ivs"	,"“", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_pound_ivs"	,"£", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_dollar_ivs"	,"$", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_hat_ivs"		,"^", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_amper_ivs"	,"&", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_obrace_ivs"	,"{", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_cbrace_ivs"	,"}", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_semi_ivs"	,";", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_pipe_ivs"	,"|", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_slash_ivs"	,"\\", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_osb_ivs"		,"[", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_csb_ivs"		,"]", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_eq_ivs"		,"=", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_lt_ivs"		,"<", CPXAxisTitle) 
CPXAxisTitle<-sub("ivs_gt_ivs"		,">", CPXAxisTitle)  
CPXAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", CPXAxisTitle) 


if(covariateModel != "NULL")
{

# Additional characters included Aug 2010 (STB)
XAxisTitle<-sub("ivs_tilde_ivs"	,"~", XAxisTitle) 
XAxisTitle<-sub("ivs_star_ivs"	,"*", XAxisTitle) 
XAxisTitle<-sub("ivs_plus_ivs"	,"+", XAxisTitle) 

XAxisTitle<-sub("ivs_sp_ivs"	," ", XAxisTitle) 
XAxisTitle<-sub("ivs_ob_ivs"	,"(", XAxisTitle) 
XAxisTitle<-sub("ivs_cb_ivs"	,")", XAxisTitle) 
XAxisTitle<-sub("ivs_div_ivs"	,"/", XAxisTitle) 
XAxisTitle<-sub("ivs_pc_ivs"	,"%", XAxisTitle) 
XAxisTitle<-sub("ivs_hash_ivs"	,"#", XAxisTitle) 
XAxisTitle<-sub("ivs_pt_ivs"	,".", XAxisTitle) 
XAxisTitle<-sub("ivs_hyphen_ivs","-", XAxisTitle) 
XAxisTitle<-sub("ivs_at_ivs"	,"@", XAxisTitle) 
XAxisTitle<-sub("ivs_colon_ivs"	,":", XAxisTitle) 
XAxisTitle<-sub("ivs_exclam_ivs","!", XAxisTitle) 
XAxisTitle<-sub("ivs_quote_ivs"	,"`", XAxisTitle) 
XAxisTitle<-sub("ivs_pound_ivs"	,"£", XAxisTitle) 
XAxisTitle<-sub("ivs_dollar_ivs","$", XAxisTitle) 
XAxisTitle<-sub("ivs_hat_ivs"	,"^", XAxisTitle) 
XAxisTitle<-sub("ivs_amper_ivs"	,"&", XAxisTitle) 
XAxisTitle<-sub("ivs_obrace_ivs","{", XAxisTitle) 
XAxisTitle<-sub("ivs_cbrace_ivs","}", XAxisTitle) 
XAxisTitle<-sub("ivs_semi_ivs"	,";", XAxisTitle) 
XAxisTitle<-sub("ivs_pipe_ivs"	,"|", XAxisTitle) 
XAxisTitle<-sub("ivs_slash_ivs"	,"\\", XAxisTitle) 
XAxisTitle<-sub("ivs_osb_ivs"	,"[", XAxisTitle) 
XAxisTitle<-sub("ivs_csb_ivs"	,"]", XAxisTitle) 
XAxisTitle<-sub("ivs_eq_ivs"	,"=", XAxisTitle) 
XAxisTitle<-sub("ivs_lt_ivs"	,"<", XAxisTitle) 
XAxisTitle<-sub("ivs_gt_ivs"	,">", XAxisTitle)  
XAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", XAxisTitle) 
}

}




#calculating number of block and treatment factors
tempblockChanges <-strsplit(blockFactors, ",")
txtexpectedblockChanges <- c("")
for(i in 1:length(tempblockChanges[[1]])) 
{
	txtexpectedblockChanges [length(txtexpectedblockChanges )+1]=(tempblockChanges[[1]][i]) 
}
noblockfactors<-length(txtexpectedblockChanges)-1

#tempChanges <-strsplit(treatFactors, ",")
#txtexpectedChanges <- c("")
#for(i in 1:length(tempChanges[[1]])) 
#{ 
#	txtexpectedChanges [length(txtexpectedChanges )+1]=(tempChanges[[1]][i]) 
#}
#notreatfactors<-length(txtexpectedChanges)-1

#need to add this code to the single measure module code?
if (tempblockChanges == "NULL")
{
	noblockfactors = 0
}

# Testing the factorial combinations
#intindex<-length(unique(statdata$betweenwithin))
#timeindex<-length(unique(eval(parse(text = paste("statdata$", timeFactor)))))
#
#ind<-1
#for (i in 1:notreatfactors)
#{
#	ind=ind*length(unique(eval(parse(text = paste("statdata$", txtexpectedChanges[i+1])))))
#}
#ind=ind*timeindex
#
#if(intindex != ind)
#{
#	add<-c("Unfortunately not all combinations of the levels of the treatment factors are present, or not all combinations are present at each level of the repeated factor, in the experimental design. We recommend you manually create a new factor corresponding to the combinations of the levels of the treatment factors or remove the 
#incomplete levels of the repeated factors from the analysis.")
#	HTML.title("</bf> ", HR=2, align="left")
#	HTML.title(add, HR=0, align="left")
#	print(add)
#	quit()
#}


# Code to create varibale to test if the highest order interaction is selected

#testeffects = noblockfactors
#if(primFactor != "NULL")
#{
#testeffects = noblockfactors+1
#}
#emodel <-strsplit(effectModel2, "+", fixed = TRUE)

#emodelChanges <- c("")
#for(i in 1:length(emodel[[1]])) 
#{ 
#	emodelChanges [length(emodelChanges )+1]=(emodel[[1]][i]) 
#}
#noeffects<-length(emodelChanges)-2




#Warning
title<-c("Warning")
HTML.title(title, HR=2, align="left")

HTML.title("</bf> ", HR=2, align="left")
HTML.title("Warning: This module is currently under construction, care should be taken when considering the results. The results have not been verified.", HR=0, align="left")




#Response
title<-c("Response")
if(primFactor != "NULL")
{
	title<-paste(title, ", covariate", sep="")
}

if (dimfact >2)
{
	title<-paste(title, " and covariance structure", sep="")
}


HTML.title(title, HR=2, align="left")
add<-paste(c("The  "), resp, sep="")
add<-paste(add, " response is currently being analysed by the Within-Subject Parametric Analysis module", sep="")
if(primFactor != "NULL")
{
	add<-paste(add, c(", with  "), sep="")
	add<-paste(add, unlist(strsplit(covariateModel, "~"))[2], sep="")
	add<-paste(add, " fitted as a covariate.", sep="")
} else {
	add<-paste(add, ".", sep="")
	}

HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")

if (responseTransform != "None" | covariateTransform != "None")
{
	HTML.title("<bf> ", HR=2, align="left")
}
if (responseTransform != "None")
{
	add2<-paste(c("The response has been "), responseTransform, sep="")
	add2<-paste(add2, " transformed prior to analysis.", sep="")
	HTML.title(add2, HR=0, align="left")
}
if (covariateTransform != "None")
{
	add3<-paste(c("The covariate has been "), covariateTransform, sep="")
	add3<-paste(add3, " transformed prior to analysis.", sep="")
	HTML.title(add3, HR=0, align="left")
}





if(covariance=="Compound Symmetric" && dimfact >2)
{
	HTML.title("</bf> ", HR=2, align="left")
	add4<-c("The Linear Mixed Model analysis is using the compound symmetric covariance structure to model the within subject correlations. When using this structure you are assuming sphericity and also that the variability of responses is the same at each level of ")
	add4<-paste(add4, timeFactor, sep= "")
	add4<-paste(add4, ", see Pinherio and Bates (2002). These assumptions may not be valid in practice.", sep ="")
	HTML.title(add4, HR=0, align="left")
}
if(covariance=="Autoregressive(1)" && dimfact >2)
{
	HTML.title("</bf> ", HR=2, align="left")
	add4<-c("The Linear Mixed Model analysis is using the first order autoregressive covariance structure to model the within subject correlations. When using this structure you are assuming the levels of ")
	add4<-paste(add4, timeFactor, sep= "")
	add4<-paste(add4, " are equally spaced and also that the variability of responses are the same at each level of ", sep= "")
	add4<-paste(add4, timeFactor, sep= "")
	add4<-paste(add4, ", see Pinherio and Bates (2002). These assumptions may not be valid in practice.", sep= "")
	HTML.title(add4, HR=0, align="left")
	HTML.title("</bf> ", HR=2, align="left")
	add4<-c("Warning: Make sure that the levels of the repeated factor occur in the correct order in the predicted means table. If they do not then this analysis may not be valid. The autoregressive covariance structure assumes that the order of the repeated factor levels is as defined in the predicted means table.")
	HTML.title(add4, HR=0, align="left")
}
if(covariance=="Unstructured" && dimfact >2)
{
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title("<bf>The Linear Mixed Model analysis is using the unstructured covariance structure to model the within subject correlations. 
When using this structure you are estimating many parameters. If the numbers of subjects used is small then these estimates may be unreliable, see Pinherio and Bates (2002).", HR=0, align="left")
}




#if(covariance=="Compound Symmetric" && dimfact ==2)
#{
#	HTML.title("</bf> ", HR=2, align="left")
#	add4<-c("The paired t-test has been calculated under the assumption that the variability is the same in both groups.")
#	HTML.title(add4, HR=0, align="left")
#}
#if(covariance=="Autoregressive(1)" && dimfact == 2)
#{
#	HTML.title("</bf> ", HR=2, align="left")
#	add4<-c("The paired t-test has been calculated under the assumption that the variability is the same in both groups.")
#	HTML.title(add4, HR=0, align="left")
#}
#if(covariance=="Unstructured" && dimfact == 2)
#{
#
#	HTML.title("</bf> ", HR=2, align="left")
#	add4<-c("The paired t-test has been calculated under the assumption that the variability is different in each group.")
#	HTML.title(add4, HR=0, align="left")
#}






#Case profiles plot
HTMLbr()
title<-c("Categorised case profiles plot of the raw data")
if(responseTransform != "None")
{
	title<-paste(title, " (on the ", sep="")
	title<-paste(title, responseTransform, sep="")
	title<-paste(title, " scale)", sep="")
}

HTML.title(title, HR=2, align="left")
scatterPlot <- sub(".html", "scatterPlot.jpg", htmlFile)
jpeg(scatterPlot)

xyplot(eval(parse(text = paste("statdata$", resp)))~eval(parse(text = paste("statdata$", timeFactor))), groups=eval(parse(text = paste("statdata$", subjectFactor ))), type="o", data=statdata, strip=strip.custom(style=1, bg="grey95"), xlab=CPXAxisTitle, ylab = YAxisTitle)#, par.settings = list( superpose.line = list(col=lev, lwd = 1) ),auto.key = list(space = "right", text=(unique(levels(as.factor(statdata$catfact)))), points = FALSE, type="b", lines = T,cex=1), col=col) 
void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", scatterPlot), Align="centre")

HTML.title("</bf> ", HR=2, align="left")
HTML.title("</bf> Tip: Use this plot to identify possible outlier (subjects and individual observations).", HR=0, align="left")

#Covariate plot
if(covar != "NULL")
{
	

	title<-c("Covariate plot of the raw data")
	if(responseTransform != "None" || covariateTransform != "None")
	{
		title<-paste(title, " (on the transformed scale)", sep="")
	} 
	HTMLbr()
	HTML.title(title, HR=2, align="left")

	primFactor2<-eval(parse(text = paste("statdata$", timeFactor)))
#	primFactor2<-paste(eval(parse(text = paste("statdata$", primFactor))))
	primFactor2 <-as.factor(primFactor2)
	statdata<-cbind(statdata, primFactor2)




	covariate<-unlist(strsplit(covariateModel, "~"))[2]

	rows<-dim(statdata)[1]
	cols<-dim(statdata)[2]
	nlevels<-length(unique(statdata$primFactor2))

	extra<-matrix(data=NA, nrow=rows, ncol=nlevels)

	for (i in 1:nlevels)
	{
		for (j in 1:rows)
		{
			if (statdata$primFactor2[j] == unique(statdata$primFactor2)[i])
			{
				extra[j,i]<-eval(parse(text = paste("statdata$", resp)))[j]
			}
		}
	}

	newdata<-cbind(statdata, extra)
	catplotdata<-data.frame(newdata)


    for (k in 1:nlevels)
    {
        tempdata<-catplotdata
        tempdata2<-subset(tempdata, statdata$primFactor2 == unique(levels(as.factor(tempdata$primFactor2)))[k])
    }

    index<-c(1:nlevels)
    newnames<-c(colnames(statdata),index)
    colnames(catplotdata)<-newnames
    ncscatterplot3 <- sub(".html", "ncscatterplot3.jpg", htmlFile)
    jpeg(ncscatterplot3)

    #Adjusting y  axis to fit in legend
    maxresp<-max(eval(parse(text = paste("statdata$",resp))))
    minresp<-min(eval(parse(text = paste("statdata$",resp))))
    rangeresp<-maxresp-minresp
    maxob<-maxresp 
#    maxob<-maxresp + rangeresp*length(unique(levels(as.factor(statdata$primFactor2))))*0.075
    minob<-minresp - rangeresp*0.1

    cat<-c(as.factor(statdata$primFactor2))
    par(las=1) 
    plot(as.formula(covariateModel), data=catplotdata, col=cat,pch=cat,  ylim=c(minob,maxob), xlab=XAxisTitle, ylab = YAxisTitle)
    

	#Best fit line 
		for (k in 1:nlevels)
		{
			tempdata<-catplotdata
			tempdata2<-subset(tempdata, tempdata$primFactor2 == unique(levels(as.factor(tempdata$primFactor2)))[k])
			abline(lm(eval(parse(text = paste("tempdata2$",resp)))~eval(parse(text = paste("tempdata2$",covariate)))), col=k)
		}
    legend("topright", legend=levels(as.factor(tempdata$primFactor2)),cex=0.6,pch=c(1:nlevels), lty=1:nlevels,bg="white", col=c(1:nlevels))
	
	void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", ncscatterplot3), Align="centre")










#	index<-c(1:nlevels)
#	newnames<-c(colnames(statdata),index)
#	colnames(catplotdata)<-newnames
#

#
#	covariatePlot <- sub(".html", "covariatePlot.jpg", htmlFile)
#	jpeg(covariatePlot)
#	plot(as.formula(covariateModel), data=catplotdata, col=c(1:nlevels),pch=c(1:nlevels))
#	void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", covariatePlot), Align="centre")
#
#	for (k in 1:nlevels)
#	{
#		tempdata<-catplotdata
#		tempdata2<-subset(tempdata, tempdata$primFactor2 == unique(levels(tempdata$primFactor2))[k])
#		abline(lm(covariateModel, data=tempdata2), col=k)
#	}
#
#	legend("topright", legend=unique(statdata$primFactor2),cex=0.6,pch=c(1:nlevels), lty=1:nlevels,bg="white", col=c(1:nlevels))

	HTML.title("</bf> ", HR=2, align="left")
	HTML.title("<bf> Tip: Is it worth fitting the covariate? You should consider the following:", HR=0, align="left")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title("<bf> a) Is there a relationship between the response and the covariate?... It is only worth fitting the covariate if there is a strong positive (or negative) relationship between them. The lines on the plot should not be horizontal.", HR=0, align="left")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title("<bf> b) Is the relationship similar for all treatments?... The lines on the plot should be approximately parallel. ", HR=0, align="left")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title("<bf> c) Is the covariate influenced by the treatment?... We assume the covariate is not influenced by the treatment so there should be no separation of the treatment groups along the x-axis on the plot. ", HR=0, align="left")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title("<bf> These issues are discussed in more detail in Morris (1999).", HR=0, align="left")
}



# Mixed model analysis code

if(covariance=="Compound Symmetric")
{
	threewayfull<-lme(model, random=~1|subjectzzzzzz, data=statdata,correlation=corCompSymm(),  na.action = (na.omit), method = "REML")
}
if(covariance=="Autoregressive(1)")
{
	threewayfull<-lme(model, random=~1|subjectzzzzzz, correlation=corAR1(value=0.4, form=~as.numeric(Timezzz)|subjectzzzzzz, fixed =FALSE), data=statdata, na.action = (na.omit), method = "REML")
}
if(covariance=="Unstructured")
{
	threewayfull<-lme(model, random=~1|subjectzzzzzz, correlation= corSymm(form = ~ as.numeric(Timezzz) | subjectzzzzzz), weights=varIdent(form=~ 1 |as.numeric(Timezzz)), data=statdata, na.action = (na.omit), method = "REML")
}



#ANOVA Table
if(showANOVA=="Y")
{
	if (dimfact ==2)
	{
	HTMLbr()
	HTML.title("<bf>Paired t-test results", HR=2, align="left")
	} else {
		HTMLbr()
		HTML.title("<bf>Overall tests of effects", HR=2, align="left")
		}		
	temp<-anova(threewayfull)
	temp2<-(temp)
	col3<-format(round(temp2[3], 2), nsmall=2, scientific=FALSE)
	col4<-format(round(temp2[4], 4), nsmall=4, scientific=FALSE)
	ivsanova<-cbind(temp2[1],temp2[2],col3,col4)
	source2<-rownames(ivsanova)
	source3<-rownames(ivsanova)
	head<-c("numDF","denDF","F-value","p-value")
	colnames(ivsanova)<-head

	for (i in 1:(dim(ivsanova)[1]-1)) 
	{
		if (ivsanova[i,4]<0.001) 
		{
			ivsanova[i,4]<-0.001
			ivsanova[i,4]<- paste("<",ivsanova[i,4])
		}
	}

	#Remove intercept row
	ivsanova <- ivsanova[-c(1), ] 
	tempy<-gsub(pattern="Timezzz", replacement=timeFactor, rownames(ivsanova))
	rownames(ivsanova)<-tempy

	HTML(ivsanova, classfirstline="second", align="left")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title("<bf>Comment: The overall tests in this table are likelihood ratio tests.", HR=0, align="left")

	add<-paste(c("Conclusion"))
	inte<-1
	for(i in 1:(dim(ivsanova)[1]))
	{
		if (ivsanova[i,4]<= (1-sig))
		{
			if (inte==1)
			{
				inte<-inte+1
				add<-paste(add, ": At the ", sep="")

				add<-paste(add, 100*(1-sig), sep="")
				add<-paste(add, "% level", sep="")

				add<-paste(add, " there is a statistically significant overall difference between the levels of ", sep="")
				add<-paste(add, rownames(ivsanova)[i], sep="")
			} else {
				inte<-inte+1
				add<-paste(add, ", ", sep="")
				add<-paste(add, rownames(ivsanova)[i], sep="")
			}
		} 
	}
	if (inte==1)
	{
		if (dim(ivsanova)[1]>2)
		{
			add<-paste(add, ": There are no statistically significant overall differences, at the ", sep="")
			add<-paste(add, 100*(1-sig), sep="")
			add<-paste(add, "% level, ", sep="")
			add<-paste(add, "between the levels of any of the effects in the overall tests table", sep="")
		} else {
			add<-paste(add, ": There is no statistically significant overall difference, at the ", sep="")
			add<-paste(add, 100*(1-sig), sep="")
			add<-paste(add, "% level, ", sep="")			
			add<-paste(add, "between the levels of the treatment effect", sep="")	
		} 
	}		
	add<-paste(add, ". ", sep="")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(add, HR=0, align="left")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title("<bf> Tip: While it is a good idea to consider the overall tests in the above table, we should not rely on them when 
	deciding whether or not to make pairwise comparisons between the factor levels.", HR=0, align="left")

	# Warning message for degrees of freedom
	if (min(ivsanova[2])<5)
	{
		HTMLbr()
		HTML.title("<bf>Warning", HR=2, align="left")
		add<-c("Unfortunately one or more of the residual degrees of freedom in the above table are low (less than 5). This may make the estimation of the underlying variability, and hence the results of the statistical tests, unreliable. Care must be taken when assessing the results of this analysis.")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(add, HR=0, align="left")
		print(add)
	}
}




#Diagnostic plot titles
if(showPRPlot=="Y")
{
	HTMLbr()
	HTML.title("<bf>Diagnostic plots", HR=2, align="left")
} else {
		if(showNormPlot=="Y")
		{
			HTMLbr()
			HTML.title("<bf>Diagnostic plots", HR=2, align="left")
		}
	}

#Residual plots
if(showPRPlot=="Y")
{
	residualPlot <- sub(".html", "residualplot.jpg", htmlFile)
	jpeg(residualPlot)
	residplot<-cbind(predict(threewayfull, level=0),residuals(threewayfull, level=0, type="pearson"))
	rownames(residplot)<-c(1:(dim(residplot)[1]))
	residplot<-data.frame(residplot)
	colnames(residplot)<-c("Predicted","Studentized_residuals")
	plot(Studentized_residuals~Predicted,data=residplot, ylab="Externally studentized residuals", xlab="Predicted values",main="Predicted vs. residuals plot")
	abline(a=0, b=0)
	abline(a=2, b=0, col="red",lty=3) 
	abline(a=3, b=0, col="red")
	abline(a=-2, b=0, col="red",lty=3)
	abline(a=-3, b=0, col="red")

	abline(h=0, col="red", lty="dotted")
	void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", residualPlot), Align="centre")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Tip: On this plot look to see if the spread of the points increases as the predicted values increase. If so the response may need transforming.", HR=0, align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Tip: Any observation with a studentized residual less than -3 or greater than 3 (SD) should be investigated as a possible outlier.", HR=0, align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Comment: The residuals at level i are obtained by subtracting the fitted levels at that level from the response vector and dividing by the estimated within-group standard error.", HR=0, align="left")
	#The residuals at level i are obtained by subtracting the fitted levels at that level from the response vector and dividing by the estimated within-group standard error.
}

#Normality plots
if(showNormPlot=="Y")
{
	HTMLbr()
	normPlot <- sub(".html", "normplot.jpg", htmlFile)
	jpeg(normPlot)
	qqnorm(resid(threewayfull, level=0) ,main="Normal probability plot")
	qqline(resid(threewayfull, level=0), col="red", lty="dotted")
	void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", normPlot), Align="left")
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Tip: Check that the points lie along the dotted line. If not then the data may be non-normally distributed.", HR=0, align="left")
}

# Means and Planned comparisons on the main effects
if(showLSMeans=="Y")
{

statdata$betweenwithin<-as.factor(eval(parse(text = paste("statdata$", timeFactor))))

	#Code to calculate y-axis offset (lens) in LSMeans plot
	names<-levels(statdata$betweenwithin)
	index<-1
	for (i in 1:length(names))
	{
		temp<-names[i]
		temp<-as.character(unlist(strsplit(as.character(names[i]),"")))
		lens<-length(temp)
		if (lens>index)
		{
			index <-lens
		}
	}

	HTMLbr()
	CITitle<-paste("<bf>Plot of the predicted means with ",(sig*100),"% confidence intervals",sep="")
	HTML.title(CITitle, HR=2, align="left")

	meanPlot <- sub(".html", "meanplot.jpg", htmlFile)
	jpeg(meanPlot)

	# LS Means


	if(covar != "NULL")
	{
		covariate<-unlist(strsplit(covariateModel, "~"))[2]
		covc<-eval(parse(text = paste("statdata$", covariate)))
		statdata<-data.frame(cbind(statdata,covc))
		covc2<-"covc"
		modl<-eval(parse(text = paste(resp, "~", covc2," + betweenwithin")))
	} else	{
			modl<-eval(parse(text = paste(resp,"~betweenwithin")))
		}



	if(covariance=="Compound Symmetric")
	{
		test<-lme(modl, random=~1|subjectzzzzzz, data=statdata, correlation=corCompSymm(), na.action = (na.omit), method = "REML")
	}
	if(covariance=="Autoregressive(1)")
	{
		test<-lme(modl, random=~1|subjectzzzzzz, correlation=corAR1(form=~as.numeric(Timezzz)|subjectzzzzzz), data=statdata, na.action = (na.omit), method = "REML")
	}
	if(covariance=="Unstructured")
	{
		test<-lme(modl, random=~1|subjectzzzzzz, correlation= corSymm(form = ~ as.numeric(Timezzz) | subjectzzzzzz), weights=varIdent(form=~ 1 |as.numeric(Timezzz)), data=statdata, na.action = (na.omit), method = "REML")

	}
	le<-dim(anova(threewayfull))[1]
	denom<-anova(threewayfull)[le,2]

#denom<-36

	#setting up the model for contrast package
	if(covar != "NULL")
	{
		covariate<-unlist(strsplit(covariateModel, "~"))[2]
		avecov<-mean(eval(parse(text = paste("statdata$", covariate))))
		covc<-eval(parse(text = paste("statdata$", covariate)))
		statdata<-data.frame(cbind(statdata,covc))
		avecov<-mean(eval(parse(text = paste("statdata$", covariate))))
		test2<-contrast(test, a = list(betweenwithin = levels(statdata$betweenwithin),covc=avecov) )
	} else
		{
			test2<-contrast(test, a = list(betweenwithin = levels(statdata$betweenwithin)))
		}

	# Creatign the table (used by gplots package)
	len<-length(test2$Contrast)
	table<-matrix(nrow=len, ncol=2)
	for (i in 1:len)
	{
		table[i, 1]=test2$Contrast[i]
		table[i, 2]=test2$SE[i]
	}
	test<-data.frame(table)
	rows<-levels(statdata$betweenwithin)
	rownames(test)<-rows
	test$lower<-test$X1-qt(1-(1-sig)/2,denom)*test$X2
	test$upper<-test$X1+qt(1-(1-sig)/2,denom)*test$X2
	cols<-c("Mean", "SE", "lower", "upper")
	colnames(test)<-cols

	#Code for LS MEans plot
	telly<- data.frame(test)
	telly2<-data.frame(rownames(test))
	telly<-data.frame(telly,telly2)
	nametemp<-c("Mean", "SE", "Lower", "Upper", "Group")
	colnames(telly)<-nametemp
	tmp1   <- split(telly$Mean, telly$Group)
	meanss <- sapply(tmp1, mean)
	tmp2   <- split(telly$Lower, telly$Group)
	lowerss <- sapply(tmp2, mean)
	tmp3   <- split(telly$Upper, telly$Group)
	upperss <- sapply(tmp3, mean)

        par(mar=c((lens/2+2),5,1,2), las=1) 
	plotCI(x=meanss, uiw=1,liw=1,li=lowerss,ui=upperss, xaxt="n", xlab=" ", ylab = YAxisTitle, col="black", barcol="black")
	par(las=2)
	axis(side=1, at=1:length(rownames(test)), cex=0.7 , labels=names(tmp1))
	void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", meanPlot), Align="left")

	HTMLbr()
	CITitle2<-paste("<bf>Table of the predicted means with ",(sig*100),"% confidence intervals",sep="")
	HTML.title(CITitle2, HR=2, align="left")

	#Code to produce the table of LS Means
	test$lower<-format(round(test$lower,2),nsmall=2)
	test$upper<-format(round(test$upper,2),nsmall=2)
	test$Mean<-format(round(test$Mean,2),nsmall=2)
	header<-c(" ", " "," ", " ")
	tables<-rbind(header, test)
	tables<-data.frame(tables)
	colnames(tables)<-c("Mean", "SE", paste("Lower ",(sig*100),"% CI",sep=""), paste("Upper ",(sig*100),"% CI",sep=""))
	rownames(tables)<-c("Level", rownames(test))
	tables$SE<-NULL

	HTML(tables, classfirstline="second", align="left")



	if(noblockfactors == 1)
	{
		add<-c("Warning: The predicted means quoted are not adjusted for the blocking factor ")
		add<-paste(add, blockFactors, sep="")
		add<-paste(add, " and hence should be used for illustrative purposes only.", sep="")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(add, HR=0, align="left")
	}

	if(noblockfactors>1) # there is two or more blocking factors
	{
		add<-c("Warning: The predicted means quoted are not adjusted for the blocking factors ")
		for (i in 1:noblockfactors)
		{
			if (i<noblockfactors-1)
			{
    			add<-paste(add, txtexpectedblockChanges[i+1], sep="")
				add<-paste(add, ", ", sep="")
			} else	if (i<noblockfactors)
	 			{
    					add<-paste(add, txtexpectedblockChanges[i+1], sep="")
					add<-paste(add, " and ", sep="")
				} else if (i==noblockfactors)
	 			{
    					add<-paste(add, txtexpectedblockChanges[i+1], sep="")
				}
		}
		add<-paste(add, " and hence should be used for illustrative purposes only.", sep="")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(add, HR=0, align="left")
	}

	if(length(grep("\\*", model)) == 0 && length(grep("\\+", model)) == 0 && length(grep("\\+", model)) == 1) 
	{
		add2<-paste(c("The effect selected involves all treatment factors."), " ", sep="")
	} else	if (length(grep("\\*", model)) == 0 && length(grep("\\+", model)) == 0 && length(grep("\\+", model)) == 0) 
		{
			add2<-paste(c(" The effect selected involves the treatment factor."), " ", sep="")
		} # else	if (noeffects>testeffects) # the model has an interaction (need to include max interaction comment)
			#{
			#	add2<-paste(c("Warning: It is not advisable to draw statistical inferences about an effect if there is a significant higher-order interaction involving that effect. In the above plot and table we have assumed that certain higher order interactions are not significant and have removed them from the statistical model, see log for 		more details."), " ", sep="")
			#}
}

#All pairwise tests
statdata$mainEffect<-as.factor(eval(parse(text = paste("statdata$", timeFactor))))
len<-length(levels(statdata$mainEffect))
k<-1
Group1<-c(1)
Group2<-c(1)
Group6<-c(1)
Group7<-c(1)
Group8<-c(1)
Group9<-c(1)

for (i in 1: (len-1))
{
	for (j in (i+1):len)
	{
		Group1[k] = levels(statdata$mainEffect)[i]
		Group2[k] = levels(statdata$mainEffect)[j]
		k=k+1
	}
}

Gplen<-length(Group1)
Group4<-strsplit(Group1, "_.._")
Group5<-strsplit(Group2, "_.._")
txtGroup4Changes<-c(" ")
for(i in 1:length(Group4[[1]])) 
{
	txtGroup4Changes [length(txtGroup4Changes )+1]=(Group4[[1]][i]) 
}
Group4len<-length(txtGroup4Changes)-1

for (i in 1:Gplen)
	{
		Group8[i]=Group4[[i]][1]
		Group9[i]=Group4[[i]][Group4len]
		Group6[i]=Group5[[i]][1]
		Group7[i]=Group5[[i]][Group4len]

		if(Group4len >2)
		{
			for (j in 2:(Group4len-1))
			{
				Group8[i]=paste(Group8[i], " ", Group4[[i]][j])
				Group6[i]=paste(Group6[i], " ", Group5[[i]][j])
			}
		}
	}

Group3<-cbind(Group6, Group7, Group8, Group9)

#Calculate denomenator df
if(covar != "NULL")
{
	covariate<-unlist(strsplit(covariateModel, "~"))[2]
	covc<-eval(parse(text = paste("statdata$", covariate)))
	statdata<-data.frame(cbind(statdata,covc))
	covc2<-"covc"
	subject<-Args[6]
	subject2<-as.factor(eval(parse(text = paste("statdata$", subject))))
	statdata<-data.frame(cbind(statdata,subject2))
	subject3<-"subject2"
	dfmodl<-eval(parse(text = paste(resp, "~", covc2, " +", subject3," + mainEffect")))
	dfanova<-anova(lm(dfmodl, data=statdata))
	dendf<-(dfanova[dim(dfanova)[1],1])
} else	{
		subject<-Args[6]
		subject2<-as.factor(eval(parse(text = paste("statdata$", subject))))
		statdata<-data.frame(cbind(statdata,subject2))
		subject3<-"subject2"
		dfmodl<-eval(parse(text = paste(resp, "~", subject3," + mainEffect")))
		dfanova<-anova(lm(dfmodl, data=statdata))
		dendf<-(dfanova[dim(dfanova)[1],1])
	}

	if(dimfact > 1)
	{
		if(covariance=="Compound Symmetric")
		{	
			test<-lme(model, random=~1|subjectzzzzzz,correlation=corCompSymm(), data=statdata, na.action = (na.omit), method = "REML")
			mult<-glht(test, linfct=mcp(Timezzz="Tukey"), df=dendf)
			multci<-confint(mult, level=sig, calpha = univariate_calpha())
		}
		if(covariance=="Autoregressive(1)")
		{	
			test<-lme(model, random=~1|subjectzzzzzz, correlation=corAR1(form=~as.numeric(Timezzz)|subjectzzzzzz), data=statdata, na.action = (na.omit), method = "REML")

			mult<-glht(test, linfct=mcp(Timezzz="Tukey"), df=dendf)
			multci<-confint(mult, level=sig, calpha = univariate_calpha())
		}
		if(covariance=="Unstructured")
		{	
			test<-lme(model, random=~1|subjectzzzzzz, correlation= corSymm(form = ~ as.numeric(Timezzz) | subjectzzzzzz), weights=varIdent(form=~ 1 |as.numeric(Timezzz)), data=statdata, na.action = (na.omit), method = "REML")
			mult<-glht(test, linfct=mcp(Timezzz="Tukey"), df=dendf)
			multci<-confint(mult, level=sig, calpha = univariate_calpha())
		}	
		multp<-summary(mult, test=adjusted("none"))	
		pvals<-multp$test$pvalues
		sigma<-multp$test$sigma
		tablen<-length(unique(rownames(multci$confint)))
		tabs<-matrix(nrow=tablen, ncol=5)

		if(dimfact > 1)
		{
			for (i in 1:tablen)
			{
				tabs[i,1]=format(round(multci$confint[i], 2), nsmall=2, scientific=FALSE)
			}
			for (i in 1:tablen)
			{
				tabs[i,2]=format(round(multci$confint[i+tablen], 2), nsmall=2, scientific=FALSE)
			}
			for (i in 1:tablen)
			{
				tabs[i,3]=format(round(multci$confint[i+2*tablen], 2), nsmall=2, scientific=FALSE)
			}
			for (i in 1:tablen)
			{
				tabs[i,4]=format(round(sigma[i], 3), nsmall=3, scientific=FALSE)
			}
			for (i in 1:tablen)
			{
				tabs[i,5]=format(round(pvals[i], 3), nsmall=3, scientific=FALSE)
			}
			for (i in 1:tablen) 
			{
				if (pvals[i]<0.001) 
				{
					tabs[i,5]<-0.001
					tabs[i,5]<- paste("<",tabs[i,5])
				}
			}

	if(dimfact > 2)
		{
		add<-paste(c("All pairwise comparisons, without adjustment for multiplicity"))
		} else 
		{
		add<-paste(c("Pairwise comparison details"))
		}

		HTMLbr()
		HTML.title(add, HR=2, align="left")
		header<-c(" ", " "," ", " "," ")
		tabls<-rbind(header, tabs)
		rows<-rownames(multci$confint)
for (i in 1:100)
{
		rows<-sub("_.._"," ", rows, fixed=TRUE)
}
		rows<-sub(" - "," vs. ", rows, fixed=TRUE)
		rownames(tabls)<-c("Comparison", rows)
		lowerCI<-paste("   Lower ",(sig*100),"% CI   ",sep="")
		upperCI<-paste("   Upper ",(sig*100),"% CI   ",sep="")
		colnames(tabls)<-c("   Difference   ", lowerCI, upperCI, "   SEM   ", "   p-value   ")
		
		HTML(tabls, classfirstline="second", align="left")
		HTML.title("<bf> ", HR=2, align="left")
	if(dimfact > 2)
		{
		HTML.title("<bf>Tip: The p-values in this table are unadjusted for multiple comparisons. No options are available in this module to make multiple comparison adjustments because it is highly unlikely you would want to make all these pairwise comparisons. If you wish to apply a multiple comparison adjustment to these results then use the p-value adjustment module.", HR=0, align="left")
		}	
		add<-paste(c("Conclusion"))
		inte<-1
		for(i in 2:(dim(tabls)[1]))
			{
				if (tabls[i,5]<= (1-sig))
				{
					if (inte==1)
					{
						inte<-inte+1
						add<-paste(add, ": The following pairwise tests are statistically significantly different at the  ", sep="")
						add<-paste(add, 100*(1-sig), sep="")
						add<-paste(add, "% level: ", sep="")
						add<-paste(add, rownames(tabls)[i], sep="")
					} else {
						inte<-inte+1
						add<-paste(add, ", ", sep="")
						add<-paste(add, rownames(tabls)[i], sep="")
						}
				} 
			}





		if (inte==1)
			{
			if (tablen >1)
				{
				add<-paste(add, ": There are no statistically significant pairwise differences.", sep="")
				} else {
					add<-paste(add, ": The pairwise difference is not statistically significant.", sep="")
					}
			} else {
				add<-paste(add, ". ", sep="")
				}
			HTML.title("</bf> ", HR=2, align="left")
			HTML.title(add, HR=0, align="left")
			if(length(grep("\\*", model)) == 0 && length(grep("\\+", model)) == 0 && length(grep("\\+", model)) == 1) 
				{
				add2<-paste(c(" "), " ", sep="")
				HTML.title("<bf> ", HR=2, align="left")
				HTML.title(add2, HR=0, align="left")
				} else	if (length(grep("\\*", model)) == 0 && length(grep("\\+", model)) == 0 && length(grep("\\+", model)) == 0) 
				{
					add2<-paste(c(" "), " ", sep="")
					HTML.title("<bf> ", HR=2, align="left")
					HTML.title(add2, HR=0, align="left")
				} #else	if (noeffects>testeffects) # the model has an interaction (need to include max interaction comment)
					#{
					#add2<-paste(c("Warning: It is not advisable to draw statistical inferences about an effect if there is a significant higher-order interaction involving that effect. In the above table we have assumed that certain higher order interactions are not significant and have removed them from the statistical model, see log for more details."), " ", sep="")
					#HTML.title("<bf> ", HR=2, align="left")
					#HTML.title(add2, HR=0, align="left")
					#}
			if (tablen >1)
				{
					HTML.title("<bf> ", HR=2, align="left")
					HTML.title("<bf>Warning: As these tests are not adjusted for multiplicity there is a risk of false positive results. Only use the pairwise tests you planned to make a-priori, these are the so called Planned Comparisons, see Snedecor and Cochran (1989).", HR=0, align="left")
				}	
			}

}


#Analysis description

HTMLbr()
HTML.title("<bf>Analysis description", HR=2, align="left")

add<-c("The data were analysed using a ")

	if(dimfact > 2)
		{
		add<-paste(add, "repeated measures Linear Mixed Model approach, with treatment factor ", sep="")
		add<-paste(add, timeFactor, sep="")
		} else {
			add<-paste(add, "paired t-test, with treatment factor ", sep="")
			add<-paste(add, timeFactor, sep="")
			}

if (blockFactors != "NULL" && primFactor != "NULL") 
{
	add<-paste(add, ", ", sep="")
} else if (noblockfactors==1 && blockFactors != "NULL" && primFactor == "NULL") 
{
	add<-paste(add, " and ", sep="")
} else if (noblockfactors>1 && blockFactors != "NULL" && primFactor == "NULL") 
{
	add<-paste(add, ", ", sep="")
} 
	
if (noblockfactors==1 && blockFactors != "NULL") 
{
	add<-paste(add, blockFactors, sep="")
	add<-paste(add, " as a blocking factor", sep="")
} else {
	if(noblockfactors>1) # there is two or more blocking factors
	{
		for (i in 1:noblockfactors)
		{
			if (i<noblockfactors-1)
			{
    			add<-paste(add, txtexpectedblockChanges[i+1], sep="")
				add<-paste(add, ", ", sep="")
			} else	if (i<noblockfactors)
	 		{
    			add<-paste(add, txtexpectedblockChanges[i+1], sep="")
				add<-paste(add, " and ", sep="")
			} else if (i==noblockfactors)
	 		{
    			add<-paste(add, txtexpectedblockChanges[i+1], sep="")
			}
		}
		add<-paste(add, " as blocking factors", sep="")
	}
}

if (covar == "NULL") 
{
	add<-paste(add, ". ", sep="")
} else if(primFactor != "NULL")	{
	add<-paste(add, " and  ", sep="")
	add<-paste(add, unlist(strsplit(covariateModel, "~"))[2], sep="")
	add<-paste(add, " as the covariate. ", sep="")
}

if (dimfact > 2)
{
	add<-paste(add, "This was followed by Planned Comparisons on the predicted means to compare the levels of the ", sep="")
	add<-paste(add, timeFactor, sep="")
	add<-paste(add, ". ", sep="")
}
if (responseTransform != "None")
{
	add<-paste(add, " The response was ", sep="")
	add<-paste(add, responseTransform, sep="")
	add<-paste(add, " transformed prior to analysis to stabilise the variance.", sep="") 
}

if (covariateTransform != "None" && responseTransform != "None")
{
	add<-paste(add, " The covariate was also ", sep="")
	add<-paste(add, covariateTransform, sep="")
	add<-paste(add, " transformed. ", sep="")
}

if (responseTransform == "None" && covariateTransform != "None")
{
	add<-paste(add, " The covariate was ", sep="")
	add<-paste(add, covariateTransform , sep="")
	add<-paste(add, " transformed prior to analysis.", sep="") 
}



HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")

if(dimfact > 2)
{
	if(covariance=="Compound Symmetric")
	{
		add2<-c("The compound symmetric covariance structure is used to model the within subject correlations. This makes the assumption that the variability of the responses is the same at each level of ")
		add2<-paste(add2, timeFactor, sep="")
	  	add2<-paste(add2, " and the correlation between responses from any pair of levels of ", sep="")
		add2<-paste(add2, timeFactor, sep="")
		add2<-paste(add2, "  is the same.")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(add2, HR=0, align="left")
	}
	if(covariance=="Autoregressive(1)")
	{
		add2<-c("The first order autoregressive covariance structure is used to model the within subject correlations. This makes the assumption that the variability of the responses is the same at each level of ")
		add2<-paste(add2, timeFactor, sep="")
	  	add2<-paste(add2, ". This structure also assumes that the correlation between responses from any pair of levels of ", sep="")
		add2<-paste(add2, timeFactor, sep="")
		add2<-paste(add2, " is related to the distance between them.")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(add2, HR=0, align="left")
	}
	if(covariance=="Unstructured")
	{
		add2<-c("The unstructured covariance structure allows the variability of the responses to be different, depending on the level of ")
		add2<-paste(add2, timeFactor, sep="")
	  	add2<-paste(add2, ". This structure also allows the correlation between responses from any pair of levels of ", sep="")
		add2<-paste(add2, timeFactor, sep="")
		add2<-paste(add2, " to be different. While this approach is the most general it should be used with care when there are few subjects, as many parameters are required to be estimated. These estimates may not be very reliable.")
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(add2, HR=0, align="left")
	}
} #else {
#
#
#	if(covariance=="Compound Symmetric")
#	{
#		add2<-c("The paired t-test has been calculated under the assumption that the variability is the same in both groups.")
#		HTML.title("</bf> ", HR=2, align="left")
#		HTML.title(add2, HR=0, align="left")
#	}
#	if(covariance=="Autoregressive(1)")
#	{
#		add2<-c("The paired t-test has been calculated under the assumption that the variability is the same in both groups.")
#		HTML.title("</bf> ", HR=2, align="left")
#		HTML.title(add2, HR=0, align="left")
#		HTML.title(add2, HR=0, align="left")
#	}
#	if(covariance=="Unstructured")
#	{
#		add2<-c("The paired t-test has been calculated under the assumption that the variability is different in the two groups.")
#		HTML.title("</bf> ", HR=2, align="left")
#		HTML.title(add2, HR=0, align="left")
#	}
#
#	}

add<-c("A full description of Linear Mixed Model theory, including information on the R nlme package used by SilveR, can be found in Venables and Ripley (2003) and Pinherio and Bates (2002).")
HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")


#References

HTMLbr()
HTML.title("<bf>Statistical references", HR=2, align="left")

if(covar != "NULL")
{
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf> Morris TR. (1999). Experimental Design and Analysis in Animal Sciences. CABI publishing. Wallingford, Oxon (UK).", HR=0, align="left")
}

HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf> Pinherio JC and Bates DM. (2000). Mixed Effects Models in S and S-Plus. Springer-Verlag. New York, Inc.", HR=0, align="left")

if (dimfact > 2)
{
	HTML.title("<bf> ", HR=2, align="left")
	HTML.title("<bf>Snedecor GW and Cochran WG. (1989). Statistical Methods. 8th edition;  Iowa State University Press, Iowa, USA.", HR=0, align="left")
}

HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf> Venables WN and Ripley BD. (2003). Modern Applied Statistics with S. 4th Edition; Springer. New York, Inc.", HR=0, align="left")



HTMLbr()
HTML.title("<bf>R references", HR=2, align="left")
	
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>   R Development Core Team (2008). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.", HR=0, align="left")

#mtvnorm
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>   Alan Genz, Frank Bretz, Torsten Hothorn with contributions by Tetsuhisa Miwa, Xuefei Mi, Friedrich Leisch and Fabian Scheipl	(2008). mvtnorm: Multivariate Normal and t Distributions. R package version 0.9-0.
	", HR=0, align="left")

#lattice
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf> 
	Deepayan Sarkar (2009). lattice: Lattice Graphics. R package version 0.17-22. http://CRAN.R-project.org/package=lattice
	", HR=0, align="left")

#gplots
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
 	Gregory R. Warnes. Includes R source code and/or documentation contributed by (in alphabetical order): Ben Bolker, Lodewijk Bonebakker, Robert Gentleman, Wolfgang Huber, Andy Liaw, Thomas Lumley, Martin Maechler, Arni Magnusson, Steffen Moeller, Marc Schwartz and Bill Venables (2009). gplots: Various R programming tools for plotting data. R package version 2.7.1. http://CRAN.R-project.org/package=gplots
	", HR=0, align="left")

#gtools
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
 	Gregory R. Warnes. Includes R source code and/or documentation contributed by Ben Bolker and Thomas Lumley (2009). gtools: Various R programming tools. R package version 2.6.1. http://CRAN.R-project.org/package=gtools
	", HR=0, align="left")

#gdata
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
 	Gregory R. Warnes, Gregor Gorjanc. Includes R source code and/or documentation contributed by Ben Bolker and Thomas Lumley. (2008). gdata: Various R programming tools for data manipulation. R package version 2.4.2.
	", HR=0, align="left")
#caTools
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
 	Jarek Tuszynski (2008). caTools: Tools: moving window statistics, GIF, Base64, ROC AUC, etc.. R package version 1.9.
	", HR=0, align="left")
#nlme
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf> 
	Jose Pinheiro, Douglas Bates, Saikat DebRoy, Deepayan Sarkar and the R Core team (2008). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-90.
	", HR=0, align="left")


#R2HTML
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
	Lecoutre, Eric (2003). The R2HTML Package. R News, Vol 3. N. 3, Vienna, Austria.
	", HR=0, align="left")

#Contrast
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>   
	Max Kuhn, contributions from Steve Weston, Jed Wing and James Forester (2009). contrast: A collection of contrast methods. R package version 0.9.
	", HR=0, align="left")

#bitops
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
	S original by Steve Dutky initial R port, extensions by Martin Maechler. revised and modified by Steve Dutky (2009). bitops: Functions for Bitwise operations. R package version 1.0-4.1.
	", HR=0, align="left")
#Survival
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
	Terry Therneau and original R port by Thomas Lumley (2009). survival: Survival analysis, including penalised likelihood.. R package version 2.35-4. http://CRAN.R-project.org/package=survival
	", HR=0, align="left")

#multcomp
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>    
	Torsten Hothorn, Frank Bretz and Peter Westfall with contributions by Richard M. Heiberger (2007). multcomp: Simultaneous Inference for	General Linear Hypotheses. R package version 0.991-8.
	", HR=0, align="left")













